﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace ODESolver
{
    public static partial class OdeSolver
    {
        /// <summary>
        /// Integrate on time T by Runge-Kutta 4  
        /// </summary>
        /// <param name="f">System of Ordinary Diff Equations</param>
        /// <param name="T">end point of integration</param>
        /// <param name="Tstep">time step</param>
        /// <param name="init">vector of initial state</param>
        /// <param name="arg">parametrs of ODESys</param>
        /// <returns></returns>
        static private double[] RungeKuttaOnTime(ODESys f, double T, double Tstep, double[] init, double[] arg)
        {
            //throw new NotImplementedException();

            int N = init.Length;
            double[] Y = new double[N];
            double t = 0;//current t
            init.CopyTo(Y, 0);

            //while (t <= (T+Tstep/2))
            while (t < T )
            {
                Y = RungeKuttaOneStep(f, t, Tstep, Y, arg);
                t += Tstep;
            }

            return Y;
        }

        static private double[,] RungeKuttaTrajectory(ODESys f, double T, double Tstep, double[] init, double[] arg)
        {
            //throw new NotImplementedException();
            int N = init.Length;
            double tmp =T /Tstep;            
            int Nsteps=(int)Math.Truncate(tmp); 
            double[,] Y = new double[Nsteps,N];
            double[] Ytemp = new double[N];
            double t = 0;//current t
            init.CopyTo(Ytemp, 0);


            for (int i = 0; i < Nsteps; i++)
            {                
                for (int j = 0; j < N; j++ )
                    Y[i, j] = Ytemp[j];

                Ytemp = RungeKuttaOneStep(f, t, Tstep, Ytemp, arg);
                t += Tstep;
            }

            return Y;
        }

        static private double[] RungeKuttaOneStep(ODESys f, double t, double dT, double[] vector, double[] arg)
        {

            #region fields
            ///<summary>
            ///Dimention of the system
            ///</summary>
            int N=vector.Length;
            ///<summary>
            ///Vector of solution
            ///</summary>
            double[] Y = new double[N];
            ///<summary>
            ///Temporary vector
            ///</summary>
            double[] Ytemp = new double[N];
            ///<summary>
            ///1st point
            ///</summary>
            double[] Y1 = new double[N];
            ///<summary>
            ///2nd point
            ///</summary>
            double[] Y2 = new double[N];
            ///<summary>
            ///3rd point
            ///</summary>
            double[] Y3 = new double[N];
            ///<summary>
            ///4th point
            ///</summary>
            double[] Y4 = new double[N];

            int i;
            #endregion

            vector.CopyTo(Y, 0);

            for (i = 0; i < N; i++)
                Y1[i] = f(i, t, Y, arg);      // рассчитать Y1

            for (i = 0; i < N; i++)
                Ytemp[i] = Y[i] + Y1[i] * (dT / 2.0);
            for (i = 0; i < N; i++)
                Y2[i] = f(i, t + dT * 0.5, Ytemp, arg);

            for (i = 0; i < N; i++)
                Ytemp[i] = Y[i] + Y2[i] * (dT / 2.0);
            for (i = 0; i < N; i++)
                Y3[i] = f(i, t + dT * 0.5, Ytemp, arg);

            for (i = 0; i < N; i++)
                Ytemp[i] = Y[i] + Y3[i] * (dT);
            for (i = 0; i < N; i++)
                Y4[i] = f(i, t + dT, Ytemp, arg);

            for (i = 0; i < N; i++)
            {
                Y[i] = Y[i] + dT / 6.0 * (Y1[i] + 2.0 * Y2[i] + 2.0 * Y3[i] + Y4[i]); // рассчитать решение на новом шаге
            }

            return Y;
        }

    }
}
